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Superoperator Many-body Theory of Molecular Currents: Non-equilibrium Green 

Functions in Real Time 
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The electric conductance of a molecular junction is calculated by recasting the Keldysh formalism 
in Liouville space. Dyson equations for nonequilibrium many body Green's functions (NEGF) are 
derived directly in real (physical) time. The various NEGFs appear naturally in the theory as time 
ordered products of superoperators, while the Keldysh forward/backward time loop is avoided. 
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I. Introduction 

Recent advances in the fabrication and measurments of nanoscale devices have lead to a considerable interest in 
nonequilibrium current carrying states of single molecules. The tunneling of electrons between two metals separated 
by a thin oxide layer was first observed experimentally by Giaever JJ and later by others . Vibrational resonances 
can be observed for molecules absorbed at the metal-oxide interface by analyzing the tunneling current as a function 
of the appHed bias0,l3|- More recent development of Scanning TunneHng Microscopy (STM) lead to a direct, real 
space determination of surface structures. A metal tip is brought near the surface so that tunneling resistence 
is measurable. A contour map of the surface is obtained by recording the tunneling resistance as the tip scans the 
surface. The tunneling electrons interact and may exchange energy with the nuclear degrees of freedom of the absorbed 
molecule. This opens up inelastic channels for electron transmission from tip to the surface leading to the ineastic 
electron tunneling (lET). lET may play an important role in manipulating molecules with STM|^, Recently, the 
lET was combined with STM for the chemical analysis of a single absorbed molecule with atomic spatial resolution 
0, 13 . Recent advances in the theory of STM are reviewed in Ref. [^. 

Electron tunneling was first analyzed by Bardeen and Cohen et a l llOLIllI using the perturbative transfer Hamilto- 
nian (TH) approach and more recently by many other authors [T^ Il3l Il4j . Although the TH gives, in most cases, a 
good description of the observed effects, it lacks a firm first principles theoretical basis and does not account properly 
for many-body effects |T5l |. An improved form of TH |l6j | that involved energy dependent transfer matrix elements, 
was used to incorporate many body effects. However this model does not describe the elecron-phonon interaction 
properlvfr^. 

A many-body non-equilibrium Green functions (NEGF) formulation of electron tunneling was proposed by Caroli 
et al [3. The NEGF theory was originated by Schwinger and Kadanoff and Baym , and developed further by 
Keldysh and Craig This formalism involves the calculation of four basic Green functions, time ordered (C^), 
anti-time ordered (G^), lesser (G*^) and greater (G^). Additional retarded (G*") and advanced (G") Green's functions 
are defined as specific combinations of these basic functions. At equilibrium suffice it to know only the retarded or 
advanced Green functions; all other Green functions simply follow from the fluctuation-dissipation theorem that 
connects the "lesser" and "greater" with the retarded Green function through the equilibrium fermi distribution 
function {fo{E)) |23j . However, for nonequilibrium measurements, where the distribution function is not known a 
priory, one needs to solve for the various NEGFs self-consistently. 

Electronic transport in molecular wires and STM currents of single molecules have received considerable attention 
O EE EE 113) EE • Electron transport through a single molecule jsE 1^112 or a chain of several atoms js^ were 
studied. From a theoretical point of view", this is very similar to the electron tunneling in semiconductor junctions 
and various theories developed for STM [SJ, |33 can directly be applied to molecular wires. The NEGF technique 
developed for tunneling currents has been used to analyze the electron conduction through a single molecule attached 
to electrodes [3 13 [sE EE EE EE EE| • The method has also been combined with density functional theory for the 
modeling of transport in molecular devices 0, E3 ■ 

In this paper we develop a nonequilibrium superoperator Green function theory '21', "2^, "43] (NESGF) of molecular 
currents |44i] . A notable advantage of working with superoperators in the higher dimensional Liouville space ,42„ ,i45| 
is that we need only consider time ordered quantities in real (physical) time; all NEGFs show up naturally without 
introducing artificial time variables. Observables can be expressed in terms of various Liouville space pathwayslhSP) 
[4^ . The ordinary (causal) response function which represents the density response to an external field is one particular 
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combination of these LSPs. Other combinations represent the spontaneous density fluctuations and the response of 
these fluctuations to the external field [isl l46l| . A simple time ordering operation of superoperators in real time 
is all it takes to derive the non-equiliubrium theory avoiding the Keldysh loop or Matsubara imaginary time. The 
NESGF theory provides new physical insights into the mechanism of the current. It can also be more naturally used 
to interpret time domain experiments involving external pulses. 

In Sec. II, we give a brief introducion to the superoperator formalism and recast the NEGF theory in terms of the 
superoperator Green functions. Starting from the microscopic definitions for various non-equilibrium superoperator 
Green functions (NESGF), we construct dynamical equations of motion and obtain the Dyson matrix equation of 
Keldysh which couples the various NESGFs. In Sec. Ill, we apply the NESGF theory to the cunduction of a molecular 
junction. In Sec. IV we end with a discussion. 



II. Dyson Equations for Superoperator Green Functions 

We consider a system of externally driven electrons and phonons described by the Hamiltonian 0, |^ [s^ , 

H = Hq + Hep + Hex (1) 

where Hq represents the non-interacting electrons and phonons, 

Ho= [ dr^\r)ho{r)^{r) + / dr(/.t(r)r!o(r),/)(r). (2) 



ho{r) = {—h?/2m)V'^ is the kinetic energy and m is the electron mass. V (V'^) represent the anihilation (creation) 
operators which satisfy the fcrmi anticommutation relations, 

^/'(r)V'^(r') + V'^(r')'/'(r) = ^(r - r') 
V't(r)i/'Hr') + V'^(r')V'^(r) = 

i/'(r)i/'(r') + ?/'(r')V'(r) = (3) 

and (f) [cf)^) are boson operators with the commutaion relations, 

0(r)0t(i.') _ 0t(i.')0(r) = ^(r-r') 
(/.^(r)0t(r') - 0t(r')(^t(i.) ^ q 

0(r)0(r') - 0(r')0(r) = (4) 

The second term in Eq. denotes the electron - phonon interaction, 

Hep = /drA(r)[(/)t(r)-f 0(r)]^t(r)^(i.). (5) 



where A(r) is the coupling strength. Finally, Hex represents the coupling to a time dependent external potantial 
C(r,t), 

Hex = / dvav,t)^\v)^j{v) (6) 



We next briefly survey some properties of Liouville space superoperators that will be useful in the following deriva- 
tions 53|- The elements of the Hilbert space N x N density matrix, p{t), wee arranged as a Liouville space vector 
(bra or ket) of length iV^. Operators of iV^ x iV^ dimension in this space are denoted as superoperators. With any 
Hilbert space operator A, we associate two superoperators Al (left) and Aji (right) defined through their action on 
another operator X as, 

AlX = AX and ArX = XA. (7) 
We further define symmetric and antisymmetric combinations of these superoperators, 

A+ - ^(^L + Ar) and A_ = {Al - Ar). (8) 
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The commutator and anticommutator operations in Hilbert space can thus be implemented with a single multiplication 
by a "-" and "+" superoperators, respectively. We further introduce the Liouville space time ordering operator T. 
This is a key ingredient for extending NEGF to superoperators: when applied to a product of superoperators it 
reorders them so that time increases from right to left. We define {A{t)) = Tr{A(i)pe^}, where peq — p{t = 0) 
represents the equilibrium density matrix of the interacting system. It is straightforward to see that for any two 
operators A and B we have, 

{TA+{t)B^it')) ^0 t'>t (9) 

{T A^{t)B ^{t')) is thus always a retarded function. This follows from the definitions ©. Since a "-" superoperator 
corresponds to a commutator in Hilbert space, this implies that for t < t' , {A+{t)B^{t')) becomes a trace of a 
commutator and must vanish, i.e., 



{TA+{t)B^{t')) = Tr{S_(t')A+(i)PeJ t<t' 
1 
2 



= T;Tr{[B{t'),A{t)p,,+p,,A{t)]} = 



Similarly, it follows that the trace of two "minus" operators always vanishes, 

(TA_(i)B_(t')) = for all t and t'. (10) 

We shall make use of Eqs. © and (|10|l for discussing the retarded and advanced Green functions in Appendix D. 
Superoperator algebra was surveyed in Ref . . 

In Liouville space the density matrix, p{t) is a vector whose time dependence is given by, 

pit)^g{t,to)pito) (11) 

with the Green function, 



g{t,to) ^Texpl^-^J^ H-(T)dT|, 



(12) 



and Ti.- is the superoperator corresponding to the Hamiltonian (Eq. Note that unlike Hilbert space, where 

time dependence of the ket and the bra is governed by forward and backward time-evolution operators respectively, in 
Liouville space one keeps track of both bra and ket simultaneously and the density matrix needs only to be propagated 
forward in time ( Eq. I^lj). 

To introduce the interaction picture in Liouville space we partition = 'Hq_ +Ti' - corresponding to the non- 
interacting and interaction Hamiltonians. With this partitinaing, Eq. (|12|l can be written as, 

g{t,to)^go{t,h)gi{t,h) (la) 

where C/o represents the time evolution with respect to Hq , 

eo(i,to) =^(i-io)exp|-^Ho_(i-to)|. (14) 
Qi{iAo) is the time evolution operator in the interaction picture. 



g/(i,io) = Texp|-^^ 7i:'_(T)dT| (15) 

and TL'_ is the interaction picture representation of 7i'_ . We shall denote superoperators in the interaction picture by 
a(~), 

i„(i) = gl{tM)Ac.{to)ga{tM (16) 
where a= -t-,- or L,R. Superoperators in the Heigenberg picture will be represented by a caret, 

A^{t) = g^{t,t„)AMQ{tM) (17) 
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By adiabatic switching of the interaction TL' - starting at = ~oo we have, 

p{t) = po-^J dTgoit,T)n'-{T)p{T) (18) 

where pq — p{—oo) is the equihbrium density matrix of the non-interacting system, 



An iterative solution of Eq. H18(l yields, 



exp{-/3Ho) 

- Tr{exp(-/3i7o)} ^''^ 



p{t) = gi{t,-^)po (20) 



which can also be obtained by applying the time evolution operator H15(l to p(io) and setting to = — oo. Using Eq. 
(|20|l . the equilibrium density matrix of the interacting system can be generated from the non- interacting one by 
switching on the interactions adiabatically starting at t = —oo. The external potential is constant in time for t < 
and is assumed to be time dependent only for i > 0. We then get, 

Peq=Gi{0,-oo)po, (21) 

This adiabatic connection formula has been shown '47] to be very useful for calculating expectaion values using the 
interaction picture. In the corresponding Gellman-law expression in Hilbcrt space |i4S | there is an extra denominator 
that takes care of the phase of the wavefunction. This is not necessery in Liouville space since the density matrix 
does not acquire such a phase. 

In the Heisenberg picture, the expectaion value of an operator A^it) is given by, 

(i„(t)) =Tr{i„(t)peg}, (22) 

where peq = p{t = 0). Using Eqs. (|16|l H17|) and (|20|l . this can be recast in the interaction picture as, 

(i,(t)> =Tr{i„(t)gj(t,-oo)po} = (i„(i)a/(t,-oo))o. (23) 

Equation H23|l is a good starting point for developing a perturbation theory around the non-interacting system. 
Through Eqs. and we also define the expectation values (...) and (...)o- While the former represents the 

trace with respect to the interacting density matrix, the latter is defined with respect to the non-interacting density 
matrix. This will be used in the following. 

Corresponding to the Hilbert space electron and phonon operators, ip, ip\ (j) and 0^ we define "left" {a—L)a.nd 
"right" {a=R) superoperators, V^q,, f/'J,, and The dynamics of a superoperator, ipa, is described by the 

generalized Liouville equation, 

_ = [n-{t)J„{t)] = n-it)^Ja.{t)-Mt)n-it)^ m 

at 

where Ti- is the superoperator corresponding to the Hamiltonian given in Eq. . A similar equation can be written 
down for the phonon superoperators. In order to evaluate the commutator appearing in the RHS of Eq. (|24|) . we 
need the commutation relations of superoperators The "left" and the "right" operators always commute. Thus 
for OL ^ (3 we have, 

[V'„(r),^^(r')] = [^l(r),V'^(r')] = [^I(r),^^(r')]=0 

[</.„(r),<^^(r')] - [0l(r),4(r')] = [0^ (r), ^^(r')] = (25) 
For fermi superoperators we have, 

V'a(r)^a(r')+V'a(r')V'a(r) = 
V't(r)V.t(r')+V4(r')V't(r) = 

V'a(r)^t(i.')+v^^(r')V-t(r) = ,5(r-r'). (26) 
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Similarly for the boson operators, 

0l(r)0l(r')-4(r')<^l(r) = 

0a(r)0„(r')-0a(r')</'a(r) = 

0„(r)0l(r')-4(r')0a(r) = «„,S(r-r'). (27) 

Here -1 for a = R and unity for a = L. 

Using the commutaion relations (|25|l . (|26|l and the identity, 

X„y„...Z„, a = L,R (28) 

we can recast 7i_ in terms of the elementary field superoperators, 

H- = Ho- + H'f + (29) 



with 



Ho- = 5] K„ /'dr(7/;t(r)/io(r)^a(r)+</.l(r)Oo(r)0a(r)) 

a=L,R •' 

Hi-" = ^ «a / drA(r)<I>„(r)V't (r)V'a(r) 

a=L,R 

= E / rfrV4(r)e(r,t)^a(r) (30) 

a=LM. '' 

We next define electron and phonon superoperator Greens functions, 

G„Mrt,rY) - -■^(r7/;„(r,t)V3^(r',i')) 

D^p{vt,v't') = -■^(r$(r,t)c|t(r',t')) (31) 

As shown in Ref. 0| (see Appendix A), Gll, Grr, Glr and Grl respectively coincide with the standard Hilbert 
space time ordered G^, antitime ordered G^, lesser G^ and greater G-^ Green functions defined on a closed time 
loop. 

Using the commutaion relations (|25|l . the Heisenberg equations of motion for superoperators i/^ait) and 3>a(t), 
a = L,R, read 

iS'ta — ' = ft.(r, i)?/'a(r, i) + A(r)$a(r,i)?/'a(r, i) 

-.HK^^^^pH = r!o(r)$„(r,i)+A(r)^t(r,i)^^(r,i), (32) 

where h(r, t) = ho{r) + ^(r, t). By taking the time derivaive of the electron Green function in Eq. (|31|l and using Eqs. 
(|32|l . we obtain the equation of motion for Gafs, 

(^h^^- K^h{v,t)^Go,p{vt,v',t') = 5„^5(x-x')--^KaA(r)(r$„(r,t)?/l„(r,t)?A;f(r',i')), (33) 
and similarly for the phonon Green function, 

(^h^^ + nM-r)^ 7^„^(ri,r',0 =5a/3<^(x-x') + -^A(r)K„(T^t(j.,t)^^(r,t)$t(r',t')) (34) 

We shall denote the space and time coordinates collectively by x = r, i; thus in Eqs. H33|l and (|34|) (5(x — x') = 
5{v-Y')6{t-t'). 

Following Keldysh, we shall rearrange the superoperator Green functions in a 2 x 2 matrix G, 

G(x,x')=f^"j"'";i (35) 

V Gi?L(x,x') Gi?ij(x,x') / ^ ^ 
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and similarly the phonon Green function matrix D with elements Dap. The corresponding Green functions of the 
noninteracting system described by the Hamiltonian ^ are denoted by and D'^, respectively. These are given by, 

G%irt,r't') = (^ih^^^K^h{r,t)^ (5„^<5(x - x') 

D%irt,r't') = (^th^^+KMr)^ 5„/3<5(x-x'). (36) 

Using our matrix notation, we can recast Eqs. (|33(l . 1)34(1 in the form of Dyson equations, 

G = GO + GOj^G 

D = DO + JJOUD. (37) 

The effect of all interactions is now included in the electron (S) and phonon (n)self energies. Exact expressions for 
the self-energies are obtained by comparing Eqs. H34|l and (|33|l with Eq. 1(3 7() . 

S„^(rt,rY) = -l«„A(r)^ f dr / dri (r$l(r, i)^„(r, i)^^,(ri, r))G^/^(rir, rY) 

J J 

n„^(rt,r'<') = ■^«aA(r)^ J dr J dri(T^t (r,t)V;„(r,t)$|,,(ri,T))i?^/^(riT,rY) (38) 

Eqs. 1(36(1 . ((37(1 and 1(38(1 are exact and constitute the non-equilibrium superoperator Green function (NESGF) 
formalism. 

In order to evaluate the self energies perturbatively. we rewrite the Green functions, Eq. ((31(1 . in the interaction 
picture, 

D^p{vt,v't') = -l(r|.(r,i)|.t(r',t')e/(t,-oo))o (39) 

where Qi{t, -co) is given by Eq. ((TB|l with = —oo. Using Eqs. (fO|l . (fTB)l and if^ . the self energies can be also 
expressed in the interaction picture as, 

S„^(ri,r'i') =-lK„A(r)^ ( f dTdri{T^i{r,t)i>^{r,t)4>l,{ruT)gi{t,~^))oG-}^{riT,r't') 

n„^(rt,r't') = '-n^X{v)Y^ J J drdri (Ti/it (r, t)Via(r, t)$^, (ri, T)e,(i, -oo))o V/J^'"!^' ^'^')- (^0) 

Equations ((4U(I together with ((39(1 constitute closed form equations for the self-energies where all the averages are 
given in the interaction picture, < ... >o, with respect to the non- interacting density matrix. By expanding Qi (Eq. 
115(1 perterbatively in 7i'_ we can obtain perturbative expansion for the self-energies. Each term in the in the expansion 
can be calculated using Wick's theorem for superoperators which is given in the Appendix E. This results in a 
perturbative series in terms of the zcroth order Green functions. 

III. The Calculation of Molecular Currents 

We have applied the NESGF to study the charge conductivity of a molecular wire attached to two perfectly 
conducting leads. In the simplest approach the leads 'a' and are treated as two free electron reservoirs. Nuclear 
motions in the molecular region are described as harmonic phonons which interact with the surrounding elctronic 
structure and the environment (secondary phonons) |24| . We first recast the general Hamiltonian, Eq. (0, in a single 
electron local basis and partition it as. 



H — Hf + Hint 



(41) 
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where Hf represents the free, non-interacting electrons and phonons and with no couphng between molecule and 
leads, 

^f=Yl ^^J^l^J + ^fc^fe'^*^- + + J2 ^r,-^^l^ra- (42) 

i,j kea,b I m 

The indices (i, j) represent the electronic basis states corresponding to the molecule, k labels the electronic states 
in the leads (a and b), I denotes primary phonons which intercts with the electrons and m denotes the secondary 
phonons which are coupled to the primary phonons and constitute a thermal bath. The applied external voltage V 
maintains a chemical potential diffrence, |Ia-^J'b = eV^, between the two leads and also modifies the single electron 
energies. In addition it provides an extra term Viipjipi which is included in the zeroth order hamiltonian, Hf, by 
modifying the single electron energies. The interaction Hamiltonian is given by, 

Hm = Y {^k^1pk^i + h.c.^ + Xu^i^l^i + (43) 

The three terms represent the molecule/lead interaction, coupling of primary phonons with the molecule and the 
interaction of primary and secndary phonons, respectively. 

The total current passing through the junction can be expressed in terms of the electron Green functions and the 
corresponding self energies. At steady state it is given by (see Appendix B, Eq. I70|l . 

- IT E / £ [^fB.i^)G'Ai^) - I]H;(-)g{«(c.)] . (44) 

ij' 

The electron Green functions G^^j^ and G^j^ correspond to the free Hamiltonian, Hf, and the self-energies ^lr and 
S_RL represent the effects of all interactions (Eq. I43|l . 

contributions coming from the electron-lead (cr) and electron-phonon (S) interactions, 

^U'^) = <^U^) + ^M^) (45) 

These are given in Eqs. (|77|l and JHEJ. The self energy expressions (|86|l and (|88|) are calculated perturbatively to 
second order in the electron-phonon coupling in terms of the zeroth order Green functions (Eq. I48|l . The simplest 
expression for current is obtained by substituting Eqs. (|48|l . (|77|1 and H86() in Eq. H44|l . This zeroth order result can 
be improved by using the renormalized Green functions obtained from the self-consistent solution of Dyson equation 

®- . 

In order to solve self-consistently for the electron Green functions that appear in the current formula, Eq. I|44l) , the 
self energy is calculated under the Born approximation by replacing the zeroth order Green functions, G^^ and £>°^ 
with the corresponding renormalized Green functions, Ga[3 and -Da/3, as is commonly done in mode-coupling theories 
|50l ISlj l . This approximation sums an infinite set of non-crossing diagrams ^.52.. .53j that appear in the perturbation 
expansion of the many body Green function, Ga(}- 

Since the electron self-energy (Eq. I86|l also depends on the phonon Green function, the phonon self-energy, H^^, 
is also required for a self-consistent solution of the electron Green functions. The phonon self-energy calculated in 
Appendix C is given by, 

= + (46) 

where Yapi^) O^'i- I78|l and A^o(a;) (Eq. I88|l represent the contributions from the phonon- phonon and the electron- 
phonon interactions, respectively. 

Computing the renormalized electron and phonon Green functions and the corresponding self-energies involves the 
self-consistent solution of the following coupled equations for the Green functions. 



Grl{oj) - G%(c.)Sfli(a;)Gii(L.)-f G%(c^)Sflfl(a;)G«i(c.) 
Gll{co) = G^jL.) + G°i(u;)Eii(L.)Gii(c.)-|-G°^(u;)Ei;^X^)GM(c^) 
GHRiij) = G"j,jiiu;)+Glj,iu^)^HLiu;)GLRii^) + Glj,iu;)i:RRiu;)GHRiij). (47) 
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Similarly the equations for the phonon Green functions are obtained by replacing Gap with Dap and Tiap with T\ai3- 

SaP 



Here Green functions corresponding to the free Hamiltonian G'^I and D'^i are given by, 



Lodij - KaEij + irj 



= ^J^B^ (48) 

where we set h— 1 and 0. Eij^Ei-Ej is the energy difference between single electron ith and jth states. Vti 
denotes the molecular phonon eigenstates. 

Once the Green functions Glr, Grl and the corresponding self-energies Sla, S/ji are obtained from the self- 
consistent solution of Eqs. H47|l together with H45|l . formula (|44|) can be used to calculate the total current through 
the molecular junction. 



Discussion 



In this paper we have developed the NESGF formalism and applied it to the computation of molecular current. 
The Liouville space time ordering operator provides an elegant way for performing calculations in real time, thus 
avoiding the artificial backward and forward time evolution required in Hilbert space (Keldysh loop) . Wick's theorem 
for superoperator is used to compute the self-enrgies perturbatively to the second order in phonon-electron coupling. 
Equations (|17|l have been derived earlier by many authors 0, |^ [s^l . Recently Galperin et al have used a 
fully self-consistent solution to study the influence of different interactions on molecular conductivity for a strong 
electron-phonon coupling. The main aim of thepresent work is to demonstrate that by doing calculations in Liouville 
space one can avoid thebackward/forward time evolution (Keldysh loop) required in Hilbert space. This originates 
from the fact that in Liouville space both ket and bra evolve forward in time. Thus one can couple the system with 
two independent fields, "left" and the "right" . This property of Liouville space can be used to construct real (physical) 
time generating functionals for the non-perturbative calculation of the self-energies. 

The pres ent model [s^ls^ls^ ignores electron-electron interactions. These may be treated using the GW technique 
|54Ll55ll56| formulated in terms of the superoperators and extended to non-equilibrium situations. All non-cquilibirium 
observables can be obtained from a single generating functional in terms of "left" and "right" operators. The retarded 
(advance) Green function that describe the forward (backward) motion of the system particle can also be calculated 
in terms of the basic Green functions, Ga/3 (see Appendix D). 

The NESGF formulation can be also recast in terms of the -I- and — (rather than L/R) superoperators which 
are more directly related to observables. This is done in Appendix D. We focused on the primary quantities that 
are represented in terms of the "left" and "right" superoperators and all other quantities are obtained as the linear 
combination of these basic operators. 
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Appendix A: Superoperator Expressions for the Keldysh Green Functions 

The standard NEGF theory formulated in terms of the four Hilbert space Green functions: time ordered (G^), 
anti-time ordered (G"^), greater (G^) and lesser fG^') [2lll23| . These are defined in the Heiscnbcrg picture as, 

G^(x,x') ^ -^(TV;(x)^t(x')) 

= ~^e{t - <')(v3(x)^t(x')) + e^t' _ t)(^t(^')^3(x)) 



G^(X,X') EE -l(fV'(x)^t(^')) 
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- -'-9{t' - t)(V'(x)^t(x')) + 0(t - i')(^t(x')v3(x)) 
G>(x,x') ^ -^(^(x)^t(^')) 

G<(x,x') ^ ■^(^t(x')^(x)) (49) 

These are known as T (T) is the Hilbert space time (anti-time) ordering operator: When appHed to a product of 
operators, it reorders them in ascending (descending) ) times from right to left. 
The four Green functions that show up naturaUy in Liouville space are defined as, 

Gll(x,x') = -^(T^i(x)V;i(x')) 

GRHi^y) = -^(r^fl,(x)v;ij(x')) 

Glk(x,x') = -^(r^i(x)^Jj(x')) 

Gfii(x,x') = -^(r^flXx)V;i(x')) (50) 

T is the LiouviUe space time ordering operator, which rearanges all superoperators in increasing order of time from 
right to left. 

To establish connection between Liouville space and Hilbert space Green functions we shall convert superoperators 
back to ordinary operators ^i^]. For Glr and Grl, we obtain, 

GLfl(x,x') = -■^Tr{TV^i(x)v31,,(x')peg} 

= -^Tr{^(x)pe,^nx')} 

= -^(^t(^')y3(x)) =G<(x,x') 

GflL(x,x') = -^Tr{r^^(x)^l(x')peg} 

= -^(^(x)^t(^')) =G>(x,x') (51) 

where peq is the fully interacting many body equilibrium density matrix. 
For Gll and Grr we have two cases, 
(i). For t > t', we get, 

Gll(x,x') = -^Tr{rviL(x)^l(x')pea 



-^Tr{^(x)^t(x')p^J = -l(^(x)^t(x')) 



Gfl,fl(x,x') = -^Tr{r^fl,(x)V^l,,(x')PeJ 



h 

i 



^Tr{pe,V^t(x')^(x)} ^ _l(^t(x')^(x)) (52) 



(ii) For the reverse case, t < t' , we get. 



Gll(x,x') EE -^Tr{rv;L(x)^l(x')peJ 



ft 

-^Tr{^t(x')y3(x)^^j ^ -l(^t(x')^(^)) 



Gflfl(x,x') = --Tr{r^^(x)V;j,,(x')PeJ 



^Tr{pe,V'(x)^nx')} - -^(^(x)^^^')) (53) 
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Combining Eqs. H52() and (|53|l we can write, 



G'll(x,x') 



--Tr{r7^i(x)v31(x')peJ 



h 



9{t - i')(V^(x)V^t(x')) - 0{t' - i)(^Hx')V'(x)) 



Gm(x,x') = 



(x,x') 

i 

h 

i 



eit - tOiV-^lxOV-lx)) - 9{t' - i)(^(x)VHx')) 



h 

= G^(x,x') (54) 
Eqs. (|51() and (|54() establish the equivalence of Hilbert and Liouville space Green functions and they can be summerized 



GLi(x,x') = G^(x,x'), G^fl(x,x') =G^(x,x') 
GLii(x,x') = G<(x,x'), Giii(x,x') =G>(x,x') 



(55) 



Appendix B: Superoperator Green Function Expression for the Current 

In this Appendix we present a formal microscopic derivation for the current flowing through a conductor. The 
conductor could be a molecule or a metal or any conducting material attached to two electrodes held at two different 
potentials. 

In Hilbert space the charge current-density is given by, 



ieh , ' 
'2m 



Vi^(r,i)VV3(r,t) - (VV;^(r,t))V;(r,t) ) 



(56) 



where e and m are the electron charge and mass, respectively. Eq. H56|l can be also expressed in a slightly modified 
form as, 



j(r,t) 



ieh r 
2m 



((V-V')V;Ur,t)^(r',t')) 



(57) 



where V represents the derivative with respect to r'. 

Using relations ()55|l the current density can be expressed in terms of the superoperator Green function as, 



j(r,t) 



2m 



[(V-V')GLfl(rt,rY: 



(58) 



At steady state, the Green functions only depend on the time difference {t — t') and the total current density (J^) 
becomes time independent. Transforming to the frequency (energy) domain, the current density per unit energy is. 



j(r, i?) = [(V - V)GLB.irv\E)l,^^ , 



and the total current density 



JT(r) 



dE 
2^ 



(59) 



(60) 



Eq. H60|) provides a recipe for calculating the current profile across the conductor once the Green function Glu is 
known from the self-consistent solution of the Dyson equation. For computing the total current passing through the 
conductor, Eq. H59() can be expressed in the form of Eq. (|44|l . In order to get the total current per unit energy 
{It{E)) passing between electrode/conductor we need to integrate the current density over the surface area of the 
conductor-electrode contact. 

(61) 



It{E) = jir,E)-hdS^ /V • j(r, E)dr 
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where n is the unit vector normal to surface S. Substituting into Eq. (|61() from (|59|l . we get 

It{E) = -I^Tr [(V2 - V^)GLR{rr',E)\ . (62) 
Zm I J 

In general, a conductor-electrode system can be described by the Hamiltonain 

H = Ho + H,^t (63) 

where Hq represents the non-interacting part, 

i/o = / dvi:\v)ho{v)4,{v) (64) 



where ha(r)= — and all the interaction terms (conductor-electrode, electron-phonon)are included in Hint- The 

total current per unit energy, Eq. (|62|l . is 

It{E) = -^Tr [(/io(r) - K{v'))Glr{vv' , E)] . (65) 

The Dyson equations for the retarded Green function (see Appendix D, Eq. I96|l . in frequency (energy) can be expressed 
in the matrix form as, 

haGr = EGr - I - (JrGr (66) 

where I is the identity matrix and is the retarded self-energy, Eq. H98(l . E=huj is a number. Henceforth we write 
all the expressins in the matrix notation. Taking the complex conjugate of (|66|l . we obtain Dyson equation for the 
advanced Green function, 

Gaho =EGa-l- Gada (67) 

with the corresponding advanced self-energy, S^. From the matrix Dyson equation H96|l . we also have the relation, 

GlR. — Gr^LRpa (68) 

Using the relations (|66|I - H68() . it is easy to see that 

hoGLR — GLR.ho — Glr'^u + Gr^LR — ^t-GlR. — ^LRpa (69) 

Substituting this in Eq. H65() . the total current per unit energy becomes, 

It{E) = '^Tr[^LR{E)GRL{E) - ^rl{E)Glr{E). (70) 

Where a factor of 2 is introduce to account for the spin degeneracy. 

We have calculated the total current in real space. In practice the Green functions and the self-energy matrices are 
calculated in an electronic basis (i,j). The total current through the conductor is obtained by integrating Eq. H7()|l 
over energy resulting in Eq. (j44(l . 



Appendix C: Self-energies for Superoperator Green Functions 

The basic quantities required for describing the coupled molecule-lead system are the one particle electron and the 
phonon Green functions. Following the steps outlined in Sec. II, the time development for various superoperators ( 
Heisenberg equations) is (all primed indices should be summed over), 

d . - „ . 

d ^ 
d - 
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where ^ma = (t^ma + 4'ma- Using Eqs. H71|) it is straightforward to write the matrix Dyson equation (|37|l for the 
electron and phonon Green functions defined as, 



G^^(x,x') = -l(r^,„(r,t)V'j^(r',i')) 
with the corresponding self-energy matrix elements, 



(72) 



- -l^o^Y. [dr 

13'.]' 



Y.yk'^{r^k'o.{t)^,fi.{T)) 



G^" (r,o-s:^(i,o+'-:i:.(M') 



nS(i,t') 



J2^l^'{T^'lJt)^^'c.m^„l,,ir)) 



D%f\T,t')^^%{t,t')+h^!i,{t,t') 



(73) 



The two terms in the electron self energy represent the contributions from the phonon-electron (S) and molecule 
- lead (fT)interactions. Similarly, the phonon self energy has contributions from the electron-phonon (A) and the 
primary-secondary phonon (7) couplings. The self energy due to the molecule-lead coupling can be calculated exactly. 
To that end we need to obtain the quantity {Tii)k' ait)'^]: piiT)) ■ By multiplying the third equation in ifTTIl by V']'/3'(''') 
from the left and from the right, taking trace and subtracting, we get (here primed indices are not summed over). 



(74) 



where gk{t) = (iftKa^ — e/c') ^ ■ Substituting expression (|74|l in Eq. (|73|l gives for the molecule-lead self energy, 

<^%{t,t') = Vk'^Vk'Jgk'{t)6{t-t') (75) 



k' ^a^b 



Similarly, the contribution to the phonon self energy from the interaction with secondary phonons can be calcualted 
exactly. 



1. 



aP = -KaSap ^ C/fm' ^^/'m'ff™' (i)'5(i - ^0 



(76) 



where g'^,{t) = {ihna§i +uJm') ^ ■ At steady state all Green functions and self-energies depend only on the time 
difference (ti — 12) and it is very convenient to express them in the frequency space. The self energy contributions due 
to molecule-lead (cr^^) and phonon-phonon (7^^) interactions, Eqs. l|75|l and {TSJ, can be represented in frequency 
space as. 



Vk'iVk'j 



■f KaUJ — tk' + «?7 



(77) 



Ulm'Uit. 



KaUJ + UJm' + iV 



(78) 
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where — > 0. However in real calculations it is a common practice to calculate self energies ctq/j and ^ap in the 
wide band approximation implying that the real parts of the self energies can be ignored and the imaginary parts are 
considered as frequency independent. Eqs. (fTTjl and l(7H|l then reduce to simpler forms, 



ii 




lal3 = 





(79) 



where T'^ = 27r^^, 14-,T4'j and T" = 2^^„, U^nVm'i'- 

The phonon contribution to the electronic self energy is obtained perturbatively in the phonon-electron coupling. 
We recast the phonon contribution (first term on the RHS of Eq. (|73|l for S^^) in the interaction picture by writing, 



(80) 



where 



C//(t, -oo) = cxp' 



X! V^'^i'a' (^)V'L' MV-^'C W 



L I' 



(81) 



Substituting (|81|l in Eq. H8U|I . expanding the exponential to first order in and using Wick's theorem for superop- 
erators Ell we obtain, 



J2 / drD'"At,r) [G°|(i,t')G°jr'(r,r+) 



(82) 



Here the superscript '0' represents the trace with respect to the non-interacting density matrix. The zeroth order 
Green functions are given in Eq. (|48() . The terms coming from the lead-molecule coupling (Vki) vanish because they 
are odd in creation and anihilation operators. Substituting (I82|l in Eq. H73|) gives for the phonon contribution to the 
self-energy. 



E%{t,t') = z^iE'*"^'-. ['*/3Ai,,i^°'/nt,t')G:i'^(t,t') 

hh 

+ 5.,6ap5{t - t') E ^o.' J drD'J^l^ {t, r)G'^X} (r, r^ 

iio.' 

In the derivation of H83|l. we have used the identity, 

/ dTY,G%{t,T)G-'f,f^{T,t') = 5^p5,,5{t~t') 



(83) 



(84) 



Similarly the contribution of the electron-phonon interaction to the phonon self energy ( second term in Eq. H73|l for 
H^g) can be obtained perturbatively. To the second order in phonon-electron coupling, we obtain. 



K'UpM = -^hY^K^np\u\l', te(t',t)G::^(t,t') +G"„^(t,<+)G°«(i',i'+) 



(85) 



To second order in electron-phonon coupling, the electronic self energy depends on both the electron and phonon 
green functions while the phonon self energy contains only the electron Green functions. 
At steady state we shift to the frequency domain and obtain. 



S^^/^M = ihY,n^n,\i,,K, f ^D^f^iu:') G°J(c.-c.') 



(86) 
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where 



pI ^ ^f^Gl^it = 0) = z / '^GtAE) (87) 



The phonon self-energy becomes, 
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Appendix D: Dyson Equations in the +/- Representation 

In this Appendix we define the retarded and advance Green's functions and the corresponding self energies and 
relate them to the basic Green functions and self energies obtained in Appendix C. From definitions Q, the Liouville 
space retarded (Gr) and advance (Ga) Green functions are defined as, 

Gl?{t,t')^-UT^j.,+ {t)ij]_{t')) (89) 



Gjf(t,i')^-^(TV^.-(OV']+(i)). 
We further introduce the correlation function, 

Gl^{t,t')^-=-{T^,+ {t')^]4t)) 



2i 



(90) 



(91) 



It follows from Eq. H1Q(I that there are only three Green functions in the +/— representation. These are given by 
Eqs. H89|) - (|91|l . Using Eq. these can be represented in terms of the basic Green functions H31|l as. 



G'iitX) 
G^{t,t') 



1 



Gtdt, t') - G'l^it, t') + Gl^{t, t') - Gt^it, t') 



= GUt,t')~G%{t,t') 

= \ [G'Lit, t') - Gl^{t, t') G%^{t, t') + Gt^{t, t') 
= ~GH«(t, t') + Gt^it, t') = Gtdt^ t') - G%,it, t') 



Gldt^ t') + Gtj^it, t') + Gtnit, t') + gH Jt, t') 



= G'l^{t,t') + G%„it,t') 



(92) 



where we have used the identity Gll + Grr ^ Glr + Grl which can be varified using Eq. (|l()|l . A Dyson equation 
corresponding to G^, Ga and Gc can be obtained from Eqs. H37|) using unitary transformation, 



G — S GS 



where G represents the matrix 



and 



G = 



Ga 
Gr Gc 



S 



1 



(93) 



(94) 



(95) 
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The transformed Dyson equation H37|l reads 

and the corresponding self energy matix reduces to 

with the matrix elements given by, 



^-Jc E^ 

E„ 



(96) 
(97) 



^RR 
^RR 



{t,t') 
{t,t') 



^'IrM 



(98) 

Similar relations also hold for the phonon Greens functions and self energies. 

Using (|98|l and (|92() , the retarded self energies for electron and phonon Green functions (retarded) coming from the 
electron-phonon coupling is obtained as, 

El^icu) = ^nY,Xl,Al, I '^^\Df{u:')Gf'{u-uj') + Df{u:')G%{u:-u:') 



w 

+ Df;,{u:')Gf\u-u:') 



duj' 



+ G'','^iu'){G°jliiu^-U^') + Ga{LU-Uj')] 



(99) 

Similarly the retarded self energies due to the lead and secondary phonons can be written in the wide band limit as 



'-r^ and 7; 



If"' 



(100) 



where F*-' includes contributions from both the leads a and b, ie, V^^— T'''J+ 



Appendix E: Wick's Theorem for Superoperators 

Wick's theorem for superoperators was formulated in Ref. Using Eqs. ijHl and H27() . it can be shown that similar 
to the L and R superoperators, the commutator of "+" and "-" boson superoperators are also numbers. Thus boson 
superoperators follow Gaussian statistics and Wick's theorem holds for both the L, R and "+", "-" representations. 
However for fermi superoperators life is more complicated. The anticommutator corresponding to only the "left" or 
the "right" fermi superoperators are numbers but that for the "left" and "right" superoperators, in general, is not a 
number. Thus the fermi superoperators are not Gaussian. However, since the left and right superoperators always 
commute, the following Wick's theorem can be applied to the time ordered product of any number of " left" and 
"right" superoperators, e.g., 

p 

Here ^i„i/„ , Vn — L, R, represents either a boson or a fermion superoperator. iaVa-^-iqi^q is a permutation of iiz^i...i„j/„ 
and sum on p runs over all possible permutations, keeping the time ordering. In case of fermions, each term should be 
multiplied by (—1)^, where P is the number of permutations of superoperators required to put them into a particular 
order. Only permutaions among either "left" or among "right" superoperators comit in P. The permutations among 
"L" and "R" operators leave the product unchanged. 
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